Hubbard model calculations of phase separation in optical lattices 
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Antiferromagnetic, Mott insulator, d-wave and gossamer superfluid phases are calculated for 2D 
square lattices from the extended Hubbard (t-J-U) model using the Gutzwiller projection method 
and renormalized mean field theory. Phase separation between antiferromagnetic and d-wave super- 
fluid phases is found near half fllling when the on-site repulsion exceeds C/;J,7.3t, and coincides with 
a first order transition in the double occupancy. Phase separation is thus predicted for 2D optical 
lattices with ultracold Fermi atoms whereas it is inhibited in cuprates by Coulomb frustration which 
instead may lead to stripes. In a confined optical lattice the resulting density distribution is discon- 
tinuous an with extended Mott plateau which enhances the antiferromagnetic phase but suppresses 
the superfluid phase. Observation of Mott insulator, antiferromagnetic, stripe and superfluid phases 
in density and momentum distributions and correlations is discussed. 

PACS numbers: 03.75.Ss, 03.75.Lm, 05.30.Fk, 74.25.Ha, 74.72.-h 



Ultra-cold atomic Fermi-gases present a new opportu- 
nity to study strongly correlated quantum many-particle 
systems and to emulate high temperature superconduc- 
tors (HTc). Optical lattices realize the Hubbard model, 
when the periodical lattice potential is strong enough so 
that only the lowest energy band is populated, and inter- 
actions, densities, temperatures, etc., can be tuned. Re- 
cent experiments with optical lattices have measured mo- 
mentum distributions and correlations, and have found 
superfluid phases of Bose and Fermi atoms 0,0, 01, Mott 
insulators (MI) and band insulators [Q, l7| . 

Long range Coulomb repulsion between electrons in 
cuprates prohibits phase separation (PS), which in turn 
may lead to stripes, whereas PS is allowed for atoms in 
optical lattices. The various competing phases can be 
calculated in the 2D t-J-U model within the Gutzwiller 
projection method and renormalized mean fleld theory 
(RMFT). This method approximates the strong correla- 
tions and generally agrees well with full variational Monte 
Carlo calculations [JQ. As will be shown below RMFT 
predicts PS at low doping between an antiferromagnetic 
(AF) Neel order and a d-wave superfluid (dSF) phase for 
sufficiently strong onsite repulsion U^7.3t. The amount 
of MI, AF and dSF phases, the density distribution and 
momentum correlations in optical lattices are quite dif- 
ferent from what would be expected from HTc cuprates. 

The t-J-U model was employed by Laughlin, Zhang 
and coworkers to study AF, HTc" and "gossamer su- 
perconductivity" in cuprates and organic superconduc- 
tors. Both the Hubbard and t-J models are included in 
the t-J-U Hamiltonian H ^ Hu + Ht+ 



or 



H = 



E 



it 



(1) 



where cua is the usual Fermi creation operator, cr (| 
, I) is the two hyperflne states (e.g. (^|, ^|) for "^^Na), 



n.ic = al^CLia the density, = J2crcr' o-la^cc'aia' and {ij) 
denotes nearest neighbours. U is the on-site repulsive 
interaction, t the nearest-neighbour hopping parameter 
and J the spin-spin or super-exchange coupHng. 

The t-J-U model allows for doubly occupied sites and 
thereby also MI transitions. As both are observed in op- 
tical lattices the t-J-U model is more useful as opposed 
to the t-J model which allows neither. For large U jt the 
t-J-U and Hubbard models reduce to the t-J model with 
spin-spin coupling J = At^ /U due to virtual hopping. At 
flnite J and U the t-J-U model is to some extent dou- 
ble counting with respect to the Hubbard model with 
J = 0. However, when RMFT is applied the virtual hop- 
ping is suppressed in the Hubbard model which justifles 
the explicit inclusion of the spin Hamiltonian as done 
in the t-J-U model. Also, optical "superlattices" provide 
additional spin-spin interaction and thus reaHze the 
unconstrained t-J-U model. 

The t-J-U model on a 2D square lattice at zero tem- 
perature was studied for various onsite coupling and den- 
sities in Refs. [1, [lH but mainly for J ~ 0.3 — 0.5f, 
which is the relevant case for cuprates and organic su- 
perconductors. In optical lattices the Hubbard model 
can be realized with any on-site repulsion by tuning near 
a Feshbach resonance. For strong onsite repulsion spin- 
exchange requires J = At^ /U and we will therefore study 
the t-J-U model with this parameter constraint. How- 
ever, keeping in mind that for more moderate U dou- 
ble counting may occur leading to an effective value for 
J. We include nearest-neighbour hopping and interac- 
tions only and therefore the phase diagram is symmetric 
around half fllling n—1, where n = N/M is the density 
or fllling fraction {N is the number of electrons or atoms 
in M lattices sites). 

In the Gutzwiller approximation the trial wave func- 
tion = ni[l — (1 — 9)fiiTfiii]\i^o) is a projection of the 
Hartree-Fock wave function tpQ. The variational param- 
eter g suppresses double occupancy and varies between 
corresponding to no doubly occupied sites {U oo) 
and 1 corresponding to no correlations {U — 0). In the 
Gutzwiller projection method the spatial correlations are 
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included only through the renormalization factors gt and 
gs defined such that (Hs) = gs{Hs)o and (Ht) = gt{Ht)o, 
where the expectation values are with respect to ip and 
ipQ respectively. The renormalization factors are calcu- 
lated by classical statistics (l^ . [l3| 



9s 



n-2d 



— 2n-|-?i_ 

for the spin term and for the hopping term 



(2) 



9t 



/gs 



— (x + d) 
2:+ 



'-±ix + d) + J^d 



(3) 



Here d is the double occupancy, n± — n/2±m the mean 
up and down occupation numbers for a magnetization m 
at each site, x = 1 — n the doping, and x± = 1 — n±. 
The double occupancy projection factor is given by g"^ = 
d(x + d)/(x+x^n^n^gs). 

The resulting energy is by Gutzwiller projection 



E = MUd + gt{Ht)o 



(4) 



The RMFT equations for the t-J-U model in the 
Gutzwiller projection method have been derived [HI in 
terms of the order parameters for d-wave pairing, hop- 
ping average and staggered magnetization (with com- 
mensurate nesting vector q = (tTjTt), in units where the 
lattice constant is unity) defined as 



A 



= ±A, 



^ = (-l)'(alT'^ii + "a"iT>o/2 , 



(5) 
(6) 
(7) 



respectively. The {+/—) in ((5| corresponds to a differ- 
ence between neighbour sites {ij) of one lattice unit in 
the {x/y) direction respectively. 

The resulting variational energy is [11] 



E/M = Ud-ilx-^Y.i^k+E, 



+ -5,J(A2+x') + 2.9,Jm2 



(8) 



Here, the AF and dSF couples four band energies Ej^ = 
)^ and 



-E^, with £,k 



A^ 



e/c = -{gtt+SgsJx/^hk, Ik = 2(cosfc:E+cosfcy) andr^^ = 
2(cosfcj, — cosky). Ad — 3gs J A/8 and Aaf — 2gsJm 
are the dSF and AF order parameters. The Lagrange 
multiplier p, differs from the chemical potential because 
the renormaHzation factors depend on density as will be 
discussed later. 

Varying the free energy F{A, x, m, d, /i) = E—jlN with 
respect to its five parameters leads to the "gap" equations 
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Figure 1: Order parameters vs. doping x = \ ~ n for J jt 
U/U = 0.1, 1/3, 0.5, 2/3 in RMFT for the 2D t-J-U model. 



for the pairing gap, hopping average and density 
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The gap equations that determines the magnetization m 
and double occupancy d [ll| are more complicated since 
the renormalization factors gs and gt also depend on m 
and d. 

The resulting order parameters for the t-J-U model 
constrained to the Hubbard or t-J model limit J — At'^/U 
are presented in Fig. 1 for U/t =40, 12, 8 and 6 corre- 
sponding to J/t = 0.1, 1/3, 0.5 and 2/3 respectively. At 
half filling (x=0) we find a first order phase transition in 
which the double occupancy jumps from zero to a finite 
value for J > 0.55t {U < 7.3t) similar to Refs. 
The ground state remains an AF for smaller U without 
any dSF. This is in agreement with expectation for the 2D 
Hubbard model but in contrast to the results for a fixed 
J = </3 where the AF vanishes for U > 7t and a dSF 
appears for [/ > 5.3t such that a coexisting AF and dSF 
at half filling (a "gossamer superfiuid") exists between 
5.3 < U/t < 7. The reason for this difference traces back 
to the larger value for J = At^ /U in our calculation which 
stabilizes the AF phase and destabilizes dSF because the 
corresponding two order parameters compete. The first 
order transition in the double occupancy is, however, ro- 
bust in the sense that it is driven by the on-site repulsion 
and remains at a critical value ?7 ~ 7 — 9t even when 
J is allowed to vary independently from the J = 4t^/J7 
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Figure 2: Phase diagram for the t-J-U model under the con- 
straint J = 44^/17. AF and PS occur for below respective 
curves. Above dotted curve the dSF gap exceeds Age > 0.2t. 
At half filling the first order transition to double occupancy 
above J > 0.55t is complementary to PS (see text). 



constraint. The case J = O.bt (U = 8t) is close to the 
critical coupling and we find a transition in d at low dop- 
ing as shown in Fig. 1 that also afi^ects the other order 
parameters m and Age- 

The gap depends sensitively on coupling, magnetiza- 
tion and double occupancy. If m = d = th e gap equa- 
tion gives A = X = Y^cos^ kx + cos^ ky/ (2-\/2M) ~ 
0.339 at half fiUing. At low densities the gap equation 
simplifies because there is no AF order Aa/ = 0. The 
gap equation depends on the pairing to quasi-particle 
energy ratio 2Ad/€k=o = VA, where V = l/(x + 
8gtt/3gsJ) — SJ/At is the effective pairing coupling di- 
vided by the bandwidth. Also the Fermi su rface be- 
comes circular with Fermi momentum kp — V2nn and 
A = i^F ~ 4)(fft + igsJx/^)- The r.h.s. of the gap equa- 
tion ^ can then be calculated analytically to leading 
logarithmic order in the effective coupling and to leading 
orders in the density. The resulting gap becomes 



1 



exp 



2(77 ^ ^" 

7m^ V 



cin — C2n 



(12) 



The higher order corrections in density can be calculated 
numerically: cq — 0.27, ci ~ 0.57 and C2 ~ 0.09. 

The AF order parameter is defined as the expec- 
tation value as in Eq. ^ but for the w.f. \'ip) in 
stead of |i/')o- It therefore also attains a renormaliza- 
tion factor niAF = \fglTn Likewise the dSF or- 

der parameter is renormalized Age = 5aA with = 

and is shown in Fig. 1. We expect that the dSF critical 
temperature is similar. 

Since the renormalization factors depend on density 
the chemical potential fi = dE/dN differs from the La- 
grange multiplier p, by 



M = A + {Ht] 



dgt 
' dn 



' dn 



Near half filling the chemical potential decreases with 
increasing density in the AF phase i.e. the energy de- 
pendence on density is concave. This is unphysical and 
signals PS to an AF phase at half filling coexisting with 
a dSC phase at a density |a;|;$,0.14 shown in Fig. 2 that 
is determined by the Maxwell construction. The PS is 
found to terminate at coupling J ~ 0.55i {U ~ 7M) 
where the double occupancy undergoes a first order tran- 
sion from zero to a finite value. Relaxing the J — At'^/U 
constraint does not change the phases or PS qualitatively 
except for gossamer dSF. 

Whereas PS is prohibited in cuprates by long range 
Coulomb repulsion between electrons, PS is permitted 
for the neutral atoms trapped in optical lattices where it 
leads to discontinuities in the density distribution n(r) 
vs. trap radius r. For a large number of trapped atoms 
in a shallow confining potential, the Thomas-Fermi ap- 
proximation applies (see, e.g. [lJ|). The total chemical 
potential is then given by the local chemical potential 
/i(n) = dE/dN and the trap potential, which is on the 
form V2r^ in most experiments. 



^l{n) + V2r^ = fi{n = 0) + V2R^ 



(14) 



(13) 



The sum must be constant over the lattice and can there- 
fore be set to its value at the edge or radius R of occu- 
pied lattice sites, which gives the r.h.s. in lfT4l) . The 
chemical potential for the dilute lattice gas in the 2D 
Hubbard model is /i(n = 0) = —At whereas just below 
half filling it Hes between ^(n = 1) = when U — Q 
and /i(n = 1) = At when U — 00. Correspondingly the 
cloud size, when PS or a MI phase appears in the cen- 
tre, lies between Rc < R < '^Rc, where Rc — 2-y/t/V2, 
and the number of trapped particles N = 2Tr n{r)rdr 
is of order Nc = ttR^. When more particles are added 
a MI plateau forms in the centre (see Fig. 3) because 
the chemical potential has a gap A^ at half filling. 
This gap and the chemical potential above half filling 
can be found from particle-hole symmetry which implies 
E{-x) = E{x) + MUx, such that /i(-x) = U-^i{x). We 
find that as U decreases so does A/i but it remains finite 
even though the double occupancy undergoes a first or- 
der transition to a finite value, i.e. the MI remains along 
with the AF phase. The approximation J = At'^/U is, 
however, only a valid for large U and the constrained t- 
J-U model may not describe the Hubbard model well for 
large J. The MI plateau remains in the trap cen tre unti l 
N^NcA^/t corresponding to a trap size R^^J A///V2. 
Increasing the number of trapped particles further en- 
forces doubly occupied sites and eventually a band insu- 
lator with n = 2 in the centre as seen in Fig. 3. 

The density distributions and MI plateaus have been 
measured experimentally for Bose atoms in optical lat- 
tices by, e.g., differentiating between singly and doubly 
occupied sites A similar technique applied to Fermi 
atoms might observe the transition in c? at C/ ~ 7.3i. 
Recently Schneider et al. [1| have measured column den- 
sities for Fermi atoms in optical lattices and find evidence 
for incompressible Mott and band insulator phases. Even 
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Figure 3: Density distributions of Fermi atoms in an optical 
lattice confined by a harmonic trap for U = 12t (J = t/3) 
filled with the number of particles N/Nc = 0.5, 1, 5, and 10. 
The density discontinuities from n ~ 1.14 and n ~ 0.86 to 
the MI at n = 1 due to PS between dSF and AF are absent 
for stripe phases (see text). 

lower temperatures are required for observing the AF and 
dSF phases and the density discontinuities due to PS. 

Spin and charge density waves in form of stripes are 
not included in the above RMFT calculations. Stripes 
are observed in several cuprates whereas numerical 
calculations are model dependent. Long range Coulomb 
frustration can explain why PS is replaced by stripes 
and an AF phase at very low doping as observed 
in cuprates. Cold atoms in optical lattices can discrim- 
inate between the PS and stripe phases since the latter 
does not have density discontinuities. 

Bragg peaks have been observed for bosons in 3D Q 
and 2D [J| optical lattices, and dips for 3D fermions in 
The Bragg peaks and dips occur in momentum cor- 
relation functions C(k, k') when the relative momentum 
is q = k — k' = Tr{nx,ny), where n^^ny are even in- 
tegers (l5l] . In an AF phase the periodicity of a given 
spin is two lattice distances and dips also appears for 



odd integers. Just as for the AF phase we can in a stripe 
phase expect charge and spin correlations as in low en- 
ergy magnetic neutron scattering p3| at incommensurate 
q — 27r(0, ±2x) and q — 7r(l, 1 ± 2x) respectively, ob- 
served as dips at these wave-numbers. However, because 
the doping x varies in the trap, the Bragg dips are dis- 
tributed over the range of values for x and are therefore 
harder to distinguish from the background. If, however, 
the four stripe periodicity occurs for l/8<x<l/4asin 
cuprates, we can expect novel Bragg dips for charge cor- 
relations at q = (0, ±7r/2) and q = (±7r/2, 0) and for spin 
correlations at q = 7r(l, 1 ± 1/4) and q = 7r(l ± 1/4, 1). 
Pairing leads to bunching between opposite momenta 
k = — k' and s-wave pairing has been observed near 
the BCS-BEC cross-over If dSF is enhanced or sup- 
pressed by stripes, the bunching due to dSF should be 
correlated correspondingly with stripe anti-bunching. 

3D optical lattices do not have a van Hove singu- 
larity at the Fermi surface at half filling as the 2D 
Hubbard model, and the 2D d-wave symmetry cx 
(cos k X cos/cy) does not generalize to 3D. Thus we do 
not expect any significant dSF but by generalizing the 
RMFT equations to 3D we find MI, AF and PS for the 
same reasons that they appear in 2D. 

In conclusion, t-J-U model RMFT calculations predict 
phase separation and density discontinuities near half fill- 
ing in 2D and 3D optical lattices when U^IM coinciding 
with a first order transition in the double occupancy at 
half filling. Observation of phase separation would indi- 
cate that long range Coulomb frustration is most likely 
the cause for spin and charge density waves and stripes in 
cuprates. Contrarily, if no PS is observed but instead a 
stripe phase near half filling in 2D and probably also 3D 
optical lattices, then Coulomb frustration is not respon- 
sible for stripes in cuprates. In either case optical lattices 
not only emulate strongly correlated systems, Hubbard 
type models and can determine the ground state phases 
such as MI, AF, PS, gossamer and dSF but can even 
determine more subtle effects from Coulomb frustration. 
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